SEISMIC  CALIBRATION  OF  GROUP  1  IMS  STATIONS  IN  EASTERN  ASIA 
FOR  IMPROVED  IDC  EVENT  LOCATION 


John  R.  Murphy  and  Jeffry  L.  Stevens 
Geophysics  Group,  Maxwell  Technologies,  Inc. 

M.  Nafi  Toksoz  and  William  L.  Rodi 
Earth  Resources  Laboratory,  MIT 

Jamil  D.  Sultanov  and  Ivan  O.  Kitov 
Institute  for  Dynamics  of  the  Geospheres 

James  F.  Lewkowicz  and  Jessie  L.  Bonner 
Weston  Geophysical  Corporation 

Sponsored  by  The  Defense  Threat  Reduction  Agency 
Arms  Control  Technology  Division 
Nuclear  Treaties  Branch 

Contract  No.  DTRA-00-C-0024 


ABSTRACT 

A  Consortium  of  institutions  which  includes  Maxwell  Technologies,  Inc.,  Massachusetts  Institute  of 
Technology’s  Earth  Resources  Laboratory  (ERL),  Weston  Geophysical  Corporation,  the  Russian  Institute  for 
Dynamics  of  the  Geospheres  (IDG),  the  Department  of  Geophysics  of  Peking  University,  and  the  Seismological 
Bureau  of  Sichuan  Province  has  been  assembled  to  conduct  a  research  program  directed  toward  the  seismic 
travel  time  calibration  of  the  30  Group  1  IMS  stations  of  eastern  Asia.  A  carefully  crafted  scientific  approach  to 
this  complex  problem  has  been  formulated  in  which  rigorous  seismological  and  statistical  analyses  of  3-D 
velocity  models  established  for  the  regions  surrounding  the  IMS  stations  of  interest  will  be  formally  integrated 
with  unique  calibration  data  sets  which  will  be  provided  by  our  international  Consortium  members  to  provide 
robust  estimates  of  source-region  specific  station  corrections  (SSSCs),  together  with  quantitative  measures  of 
the  uncertainties  in  those  estimates.  A  rationale  for  assigning  priorities  to  the  IMS  stations  to  be  investigated 
has  been  established,  and  a  staged  program  of  deliverables  has  been  formulated  in  which  preliminary  calibration 
results  for  these  priority  stations  will  be  provided  six  months  after  contract  award  and  updated  and  refined  on  a 
regular  basis  throughout  the  three  year  performance  period.  We  are  currently  assembling  a  preliminary  3-D 
velocity  model  of  the  entire  region  which  is  composed  of  a  global  background  model  on  a  5°-by-5°  grid 
(Stevens  and  McLaughlin,  1999),  supplemented  by  more  detailed  models  in  regions  where  they  are  available. 
At  the  present  time,  such  detailed  models  have  been  identified  for  a  large  portion  of  the  Former  Soviet  Union 
(FSU)  for  which  Deep  Seismic  Sounding  (DSS)  data  have  been  used  to  define  a  3-D  velocity  model  of  the  crust 
and  upper  mantle  on  a  roughly  40km-by-40km  grid  (Shchukin,  1989),  and  for  an  approximately  20°-by-20°  area 
centered  on  the  Pakistan/ Afghanistan  region  for  which  a  3-D  velocity  model  has  been  defined  on  a  l°-by-l° 
grid  (Bernard  et  al.,  1999).  Regional  phase  travel  times  through  these  3-D  velocity  models  are  being  computed 
using  the  Podvin  and  Lecomte  (1991)  finite  difference  algorithm  to  obtain  preliminary  SSSC  estimates  with 
respect  to  the  IDC  default  IASPEI91  model  for  each  of  the  30  Group  1  IMS  stations.  These  preliminary 
estimates  will  be  tested  and  refined  using  the  various  calibration  data  sets  which  are  being  assembled  for  this 
study.  These  include  a  unique  set  of  regional  arrival  time  data  at  60  FSU  permanent  network  stations  from 
some  40  Soviet  PNE  tests,  as  well  as  50  Semipalatinsk  and  Novaya  Zemlya  explosions  with  precisely  known 
locations  and  origin  times,  including  a  number  of  stations  which  are  either  at  or  very  near  to  Group  1  IMS 
station  locations.  The  resulting  SSSCs  will  then  be  thoroughly  tested  to  quantify  improvements  to  event 
locations  in  terms  of  uncertainty  and  location  errors  relative  to  the  known  locations  of  selected  regional  ground 
truth  events.  Our  ultimate  objective  is  to  provide  the  IDC  with  the  enhanced  location  capability  which  will  be 
required  to  meet  the  CTBT  monitoring  objectives  for  this  large  Group  1  IMS  station  region. 
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OBJECTIVES 


The  purpose  of  this  effort  is  to  develop  improved  3-D  velocity  models,  Site-Specific  Station  Corrections 
(SSSCs),  and  Slowness-Azimuth  Station  Corrections  (SASCs)  for  eastern  Asia,  to  demonstrate  the  effectiveness 
of  these  models  and  corrections  in  improving  locations  of  seismic  events,  to  evaluate  the  uncertainties 
associated  with  these  improved  models  and  corrections,  and  to  install  these  models  and  corrections  at  the 
Prototype  International  Data  Centre  (PIDC)  of  the  Comprehensive  Nuclear-Test-Ban  Treaty  (CTBT)  located  at 
the  Center  for  Monitoring  Research  (CMR)  and  evaluate  their  performance.  The  specific  objectives  are 
fourfold.  The  first  is  to  collect  appropriate  calibration  data  and  seismic  events  occurring  in  and  around  the 
calibration  region.  The  second  is  to  use  these  data  to  refine  the  regional  velocity  models  and  to  define  and  refine 
SSSCs  and  SASCs  for  the  30  IMS  stations  within  the  region.  The  third  is  to  use  the  new  velocity  models  and 
SSSCs  and  SASCs  with  calibration  data  [“Ground  Truth”  (GT)  data]  to  evaluate  the  location  capabilities  of  the 
system.  The  fourth  is  to  implement  the  SSSCs  and  SASCs  into  the  location  system  at  PIDC/CMR  and  to  work 
with  the  CMR  staff  to  evaluate  their  performance  in  that  operational  environment. 

RESEARCH  ACCOMPLISHMENTS 

The  determination  of  accurate  seismic  locations  for  detected  events  is  one  of  the  most  important  functions  of  the 
IDC  because  location  plays  such  a  key  role  in  CTBT  verification.  Thus,  for  example,  under  the  CTBT,  the  areas 
to  be  designated  for  potential  onsite  inspections  are  limited  to  1000  km2;  and,  consequently,  this  value  has  been 
adopted  as  the  minimum  performance  requirement  for  uncertainties  in  IDC  seismic  locations.  However,  it  has 
proven  difficult  to  demonstrate  that  this  level  of  accuracy  can  be  routinely  achieved,  particularly  for  small 
events  recorded  by  limited  numbers  of  IMS  stations.  The  principal  reason  for  this  difficulty  is  that  the  earth  is 
laterally  heterogeneous  over  a  very  broad  range  of  length  scales;  and,  therefore,  travel  time  curves  based  on 
radially  symmetric  earth  models,  such  as  the  standard  IASPEI91  model  employed  at  the  IDC,  can  be  seriously 
in  error  for  some  source  and  station  locations.  This  problem  is  particularly  acute  at  regional  distances  where  the 
propagation  paths  of  the  seismic  phases  used  in  location  are  predominantly  confined  to  the  crust  and  upper 
mantle  portions  of  the  Earth  which  exhibit  the  greatest  regional  variability.  Since  many  of  the  small  seismic 
events  which  are  of  primary  concern  in  CTBT  monitoring  will  have  to  be  located  using  regional  seismic  data, 
such  variability  places  significant  limitations  on  currently  achievable  location  accuracy.  It  follows  that  in  order 
to  establish  a  robust  CTBT  monitoring  capability,  it  will  first  be  necessary  to  calibrate  the  IMS  stations  to 
account  for  such  systematic  deviations  from  the  nominal  travel  time  curves. 

We  are  in  the  early  stages  of  a  three  year  program  which  has  the  objective  of  calibrating  the  30  Group  1  IMS 
stations  of  eastern  Asia.  The  map  locations  of  these  stations  are  shown  in  Figure  1,  where  it  is  indicated  that 
they  are  composed  of  1 1  primary  and  19  secondary  stations.  It  can  be  seen  from  Figure  1  that  regional  coverage 
circles  extending  to  2000  km  radius  surrounding  these  stations  would  encompass  most  of  the  Asian  continent, 
with  overlapping,  multiple  station  coverage  in  most  areas.  Thus,  the  area  under  investigation  is  very  large 
indeed.  Of  this  total  of  30  stations,  13  are  currently  active  IMS  stations  and  7  more  are  in  close  proximity  to 
permanent  network  stations  for  which  significant  historical  calibration  data  are  available.  This  subset  of  20 
stations  has  been  designated  as  the  “Priority  1”  stations  for  the  purposes  of  our  investigation. 

Our  overall  approach  to  this  problem  centers  on  the  formulation  of  a  3-D  velocity  model  of  the  region  which 
can  be  used  to  define  corrections  to  the  default  IASPEI91  model  currently  employed  at  the  IDC.  Thus,  an  initial 
velocity  model  is  defined  from  currently  available  information  and  is  then  refined  by  performing  joint 
tomography  and  multiple  event  locations  using  arrival  time  data  collected  for  this  region.  Once  an  optimum 
average  model  has  been  developed,  source-specific  empirical  travel  time  corrections  for  each  station  will  be 
determined  by  interpolating  between  observed  calibration  event  residuals  relative  to  this  3-D  model.  Model- 
based  travel  times  determined  by  ray  tracing  through  the  3-D  velocity  model  will  then  be  combined  with  these 
empirical  corrections  to  generate  3-D  travel  time  tables  for  each  Group  1  IMS  station,  and  these  predicted  travel 
times  will  be  differenced  against  the  corresponding  IASPEI91  travel  times  to  define  the  required  SSSCs  as  a 
function  of  source  location  around  each  station.  This  process  will  then  be  iterated  to  incorporate  data  from 
additional  calibration  events  as  they  become  available  and  the  final  models  will  be  validated  based  on  relocation 
experiments  conducted  using  travel  time  data  recorded  from  ground  truth  calibration  events  in  the  region. 


Our  preliminary  model  for  the  Group  1  region  is  a  composite  of  three  models:  the  long-period  models  of 
Stevens  and  Adams  (1999)  (model  LP),  the  preliminary  model  of  Bernard  et  al.  (1999)  for  the  Pakistan  region 
(model  OPM)  and  the  Deep  Seismic  Sounding  models  for  the  former  Soviet  Union  (model  DSS).  Each  of  these 
models  is  defined  on  a  regular  latitude/longitude  grid  of  laterally  homogeneous  cells.  The  cells  of  model  LP  are 
5  degrees  by  5  degrees  and  cover  the  entire  globe.  The  cells  of  model  OPM  are  1  degree  by  1  degree  and  cover 
24.5  degrees  N  to  43.5  degrees  N  in  latitude  and  54.5  degrees  E  to  80.5  degrees  E  in  longitude.  The  DSS  model 
cells  are  0.3333  degrees  wide  in  latitude  and  0.5  degrees  wide  in  longitude  and  cover  40  N  to  72.33  N  latitude 
and  30  to  140  degrees  E  longitude.  Many  of  the  DSS  cells,  however,  are  not  defined.  Each  cell  of  each  model  is 
associated  with  a  layered  velocity  versus  depth  profile.  The  LP  profiles  extend  to  600  km  depth.  The  OPM  and 
DSS  profiles  extend  to  30  km  below  the  Moho  depth,  which  is  generally  different  for  different  cells.  The  areas 
covered  by  these  three  models  are  identified  on  the  map  in  Figure  2,  where  the  annotations  in  the  5  by  5  degree 
background  model  denote  the  different  1-D  velocity  models  inferred  by  Stevens  and  Adams  (1999)  from  their 
tomographic  inversion  analyses  of  observed  IMS  surface  wave  data. 

The  three  models  were  assembled  in  hierarchy  as  follows.  First,  a  master  rectangular  grid  of  latitude/longitude 
cells  was  constructed  having  variable  cell  dimensions  and  containing  all  the  cell  boundaries  of  the  three  input 
models.  Each  cell  of  the  master  grid  was  assigned  a  composite  velocity  profile  defined  to  be  the  appropriate  LP 
profile  with  the  shallower  depths  replaced  by  an  OPM  profile  or  DSS  profile,  when  either  is  defined.  For  a 
small  number  of  master  cells,  both  OPM  and  DSS  profiles  are  defined;  and  we  allowed  the  latter  to  take 
precedence.  If  the  deepest  velocity  of  a  DSS  or  OPM  profile  exceeds  the  LP  velocity  at  that  depth,  the  velocities 
of  LP  layers  below  that  depth  were  raised  to  this  value.  We  point  out  that  the  master  model,  like  LP,  has  profiles 
which  vary  from  cell  to  cell  in  the  number  of  total  layers,  number  of  crustal  layers,  etc. 

The  raytracing  algorithm  requires  a  grid  of  constant-velocity,  cubic  cells  in  a  Cartesian  coordinate  system.  We 
constructed  such  a  grid  for  each  station,  extending  2000  km  north,  south,  east  and  west  of  the  station,  and  to  a 
depth  of  600  km.  The  cell  size  was  5  km.  For  each  vertical  column  of  cells  in  the  Cartesian  raytracing  grid,  we 
mapped  the  horizontal  coordinates  of  the  center  of  the  column  to  a  latitude  and  longitude,  based  on  a 
sterographic  projection  centered  at  the  station.  This  latitude  and  longitude  is  then  located  in  the  composite 
velocity  model  and  the  velocity  profile  at  this  point  is  extracted.  The  velocity  profile  is  then  translated  into  a 
uniform  5  km  layering  by  calculating  layer  velocities  that  preserve  the  vertical  travel  time  to  each  interface. 
This  translation  is  done  in  conjunction  with  an  earth-flattening  transformation,  which  extends  the  legitimacy  of 
Cartesian  raytracing  to  greater  distances  from  the  station.  The  re -layered  velocity  profile  defines  the  cell 
velocities  for  a  column  of  the  raytracing  grid. 

Ray  tracing  through  this  complex  3-D  velocity  model  is  being  carried  out  using  a  code  based  on  the  finite 
difference  algorithm  of  Podvin  and  Lecomte  (1991).  There  are  several  advantages  to  using  this  approach  in  the 
current  application.  First,  it  permits  the  inclusion  of  secondary  arrivals  in  the  location  problem;  and  secondly,  it 
improves  on  previous  modeling  methods  by  being  able  to  accommodate  very  large  velocity  contrasts  of 
arbitrary  shape,  as  high  as  10:1.  It  is  particularly  convenient  for  estimation  of  SSSCs  in  that  transmitted  and 
diffracted  body  waves,  as  well  as  head  waves,  can  be  readily  modeled  for  each  possible  hypocentral  location  in 
a  Cartesian  grid  centered  on  a  particular  station.  A  first  arrival  criterion  is  then  used  to  choose  among  the 
various  types  of  waves  to  create  the  travel  time  tables.  This  algorithm  has  been  extensively  tested  against  a 
variety  of  models,  including  that  corresponding  to  IASPEI91,  and  found  to  be  accurate  to  better  than  0.5 
seconds  over  the  entire  range  of  source  to  station  distances  of  interest  in  this  project. 

It  is  well  known  that  there  are  extensive  portions  of  the  study  area  shown  in  Figure  1  which  are  essentially 
aseismic  and,  therefore,  cannot  be  effectively  calibrated  using  earthquake  data  alone.  Fortunately,  however,  the 
U.S.S.R.  did  conduct  an  extensive  series  of  PNE  tests  at  widely  dispersed  locations  throughout  the  territories  of 
the  FSU  for  which  precise  locations  and  origin  times  are  now  available  (Sultanov  et  al.,  1999).  Travel  time  data 
from  these  explosions  were  observed  at  numerous  stations  of  the  U.S.S.R.  permanent  seismic  network,  and 
these  data  obviously  constitute  a  primary  resource  for  the  calibration  of  Group  1  IMS  stations  in  the  FSU. 
Figure  3  shows  the  map  locations  of  the  approximately  60  U.S.S.R.  permanent  network  stations  for  which  travel 
time  data  are  currently  available  from  Soviet  nuclear  tests.  Although  all  of  these  stations  can  provide  ground 
truth  data  useful  for  calibrating  subsurface  velocity  models  for  the  region,  there  is  a  subset  of  8  of  them  which 
are  either  at  or  very  near  to  Group  1  IMS  station  locations;  and  these  will  obviously  provide  very  important  data 


calibration  database  of  regional  ground  truth  data  is  available  from  over  100  nuclear  explosions.  This  fact  is 
illustrated  in  Figure  4  which  shows  the  locations  of  underground  nuclear  tests  within  about  20°  of  the 
Borovoye  station,  together  with  the  observed  P  wave  travel  time  residuals  with  respect  to  IASPEI  at  BRVK 
for  each  test  location.  It  can  be  seen  from  this  figure  that  these  regional  calibration  events  are  very  well 
distributed  about  the  BRVK  station  and  that  the  observed  travel  time  residuals  with  respect  to  IASPEI  show 
a  very  complex  pattern  varying  from  about  -8.0  to  +1.0  seconds.  It  is  probably  fair  to  say  that  there  is  no 
other  IMS  station  location  for  which  a  comparable  amount  of  well  distributed  ground  truth  travel  time  data 
is  currently  available.  Consequently,  we  have  chosen  to  use  BRVK  as  our  benchmark  station  for  the  testing 
and  evaluation  of  all  of  our  calibration  procedures. 

As  an  initial  test  of  our  modeling  procedures,  we  have  estimated  a  preliminary  SSSC  for  station  BRVK  by 
ray  tracing  through  a  3-D  DSS  velocity  model  covering  the  region  out  to  20°  around  that  station.  Since  DSS 
Pn  velocities  are  not  currently  available  for  this  entire  area,  we  have  adopted  an  upper  mantle  Pn  velocity 
model  which  is  consistent  with  the  average  observed  PNE  arrival  times  as  a  function  of  distance  at  BRVK 
for  this  initial  application.  A  grayscale  plot  of  the  resulting  surface  focus  SSSC  is  shown  in  Figure  5  where 
it  can  be  seen  that  the  preliminary  3-D  velocity  model  predicts  a  complex  pattern  of  corrections  with 
respect  to  IASPEI,  with  a  total  range  of  variation  extending  from  about  -7.9  to  +2.7  seconds.  The  extent  to 
which  this  model  predicts  the  observed  PNE  ground  truth  residuals  of  Figure  4  is  illustrated  in  Figure  6, 
which  shows  the  corresponding  residuals  computed  with  respect  to  the  travel  times  predicted  by  the  3-D 
DSS  velocity  model.  It  can  be  seen  that  the  residuals  with  respect  to  this  latter  model  are  greatly  reduced 
relative  to  those  of  Figure  4,  with  the  majority  of  residual  values  falling  in  the  range  of  about  +1.5  seconds. 
This  fact  is  illustrated  more  clearly  in  Figure  7,  where  the  travel  time  residuals  from  Figures  4  and  6  are 
compared  on  an  event  by  event  basis.  It  can  be  seen  that,  while  the  residuals  with  respect  to  the  original 
IASPEI  model  show  an  average  bias  of  approximately  -3.7  sec,  those  computed  with  respect  to  the  DSS 
model  are  essentially  unbiased,  showing  an  average  value  of  only  -0.1  seconds.  Moreover,  the  standard 
deviation  of  the  residuals  decreases  from  about  2.2  to  1.7  seconds  in  going  from  the  IASPEI  to  the  DSS 
model,  indicating  that  the  DSS  model  not  only  eliminates  the  average  bias,  but  also  accounts  for  more  of 
the  site  specific  variations  between  individual  events.  Thus,  although  further  refinements  are  clearly 
required,  the  initial  results  obtained  using  the  3-D  DSS  velocity  model  are  encouraging.  We  are  currently 
proceeding  to  apply  these  same  modeling  procedures  to  the  estimation  of  preliminary  SSSCs  for  all  20  of 
our  Priority  1  stations. 


CONCLUSIONS  AND  RECOMMENDATIONS 

We  are  in  the  early  stages  of  a  comprehensive  program  to  calibrate  the  30  Group  1  IMS  stations  of  eastern 
Asia  in  an  attempt  to  improve  IDC  seismic  location  capability  in  this  area.  A  preliminary,  composite  3-D 
velocity  model  for  the  entire  region  has  now  been  formulated  and  an  efficient  finite  difference  ray  tracing 
code  for  computing  travel  times  through  such  models  has  been  implemented  and  tested.  A  preliminary 
ground  truth  database  consisting  of  arrival  time  observations  at  approximately  60  stations  of  the  U.S.S.R. 
permanent  seismic  network  from  Soviet  PNE  tests  has  been  assembled,  including  data  at  a  number  of 
stations  which  are  either  at  or  very  near  to  Group  1  IMS  station  locations.  Among  these  is  the  Borovoye 
station  BRVK  for  which  an  exceptionally  well  distributed  ground  truth  data  set  of  regional  arrival  time 
observations  from  Soviet  PNE  events  is  available  for  testing  of  proposed  calibration  procedures.  Initial  tests 
of  a  3-D,  DSS-based  velocity  model  for  the  region  around  station  BRVK  have  been  encouraging  in  that 
they  indicate  that  such  a  model  can  account  for  virtually  all  of  the  observed  travel  time  bias  and  a 
significant  amount  of  the  site  to  site  variability  associated  with  the  predictions  of  the  IDC  default  IASPEI91 
model.  Current  efforts  are  directed  toward  the  estimation  of  preliminary  SSSCs  for  the  subset  of  20  Group 
1  IMS  stations  which  have  been  designated  as  the  highest  priority  for  calibration. 


Figure  3.  Map  locations  of  U.S.S.R.  permanent  network  stations  for  which  travel  time  data  are 
currently  available  from  Soviet  nuclear  tests. 


Figure  4.Map  locations  of  regional  distance  underground  nuclear  tests  recorded  at  IMS  station 
location  BRVK,  together  with  associated  observed  travel  time  residuals  with  respect  to  IASPEI. 
The  oval  line  marks  a  constant  epicentral  distance  of  20°  from  BRVK. 
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Figure  5.  Surface  focus  SSSC  for  the  Borovoye  station  computed  using  the  DSS  based  3-D  velocity 
model  for  the  surrounding  region.  Differential  arrival  times  are  computed  with  respect  to  the  IDC 
default  IASPEI91  model. 


-as  *+ 
sa-as  s' 


s  sj?  Is  '  WA’PXT' 


Figure  6.  P  wave  travel  time  residuals  for  PNE  ground  truth  events  recorded  at  the  BRVK  station, 
computed  with  respect  to  the  travel  times  predicted  by  the  DSS  based  3-D  velocity  model  for  the 
surrounding  region. 
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Figure  7.  Comparison  of  P  wave  travel  time  residuals  for  PNE  events  recorded  at  the  BRVK  station, 
computed  with  respect  to  the  IASPEI  91  and  3-D  DSS  velocity  models,  respectively. 
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